<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Handling functions with large features near an endpoint with tanh-sinh quadrature</title>
<link rel="stylesheet" href="../../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../../index.html" title="Math Toolkit 3.3.0">
<link rel="up" href="../double_exponential.html" title="Double-exponential quadrature">
<link rel="prev" href="de_tanh_sinh.html" title="tanh_sinh">
<link rel="next" href="de_sinh_sinh.html" title="sinh_sinh">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="de_tanh_sinh.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../double_exponential.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="de_sinh_sinh.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h3 class="title">
<a name="math_toolkit.double_exponential.de_tanh_sinh_2_arg"></a><a class="link" href="de_tanh_sinh_2_arg.html" title="Handling functions with large features near an endpoint with tanh-sinh quadrature">Handling
      functions with large features near an endpoint with tanh-sinh quadrature</a>
</h3></div></div></div>
<p>
        Tanh-sinh quadrature has a unique feature which makes it well suited to handling
        integrals with either singularities or large features of interest near one
        or both endpoints, it turns out that when we calculate and store the abscissa
        values at which we will be evaluating the function, we can equally well calculate
        the difference between the abscissa value and the nearest endpoint. This
        makes it possible to perform quadrature arbitrarily close to an endpoint,
        without suffering cancellation error. Note however, that we never actually
        reach the endpoint, so any endpoint singularity will always be excluded from
        the quadrature.
      </p>
<p>
        The tanh_sinh integration routine will use this additional information internally
        when performing range transformation, so that for example, if one end of
        the range is zero (or infinite) then our transformations will get arbitrarily
        close to the endpoint without precision loss.
      </p>
<p>
        However, there are some integrals which may have all of their area near
        <span class="emphasis"><em>both</em></span> endpoints, or else near the non-zero endpoint,
        and in that situation there is a very real risk of loss of precision. For
        example:
      </p>
<pre class="programlisting"><span class="identifier">tanh_sinh</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span> <span class="identifier">integrator</span><span class="special">;</span>
<span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span> <span class="special">/</span> <span class="special">(</span><span class="number">1</span> <span class="special">-</span> <span class="identifier">x</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">);</span> <span class="special">};</span>
<span class="keyword">double</span> <span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="number">1.0</span><span class="special">);</span>
</pre>
<p>
        Results in very low accuracy as all the area of the integral is near 1, and
        the <code class="computeroutput"><span class="number">1</span> <span class="special">-</span>
        <span class="identifier">x</span> <span class="special">*</span> <span class="identifier">x</span></code> term suffers from cancellation error
        here.
      </p>
<p>
        However, both of tanh_sinh's integration routines will automatically handle
        2 argument functors: in this case the first argument is the abscissa value
        as before, while the second is the distance to the nearest endpoint, ie
        <code class="computeroutput"><span class="identifier">a</span> <span class="special">-</span>
        <span class="identifier">x</span></code> or <code class="computeroutput"><span class="identifier">b</span>
        <span class="special">-</span> <span class="identifier">x</span></code>
        if we are integrating over (a,b). You can always differentiate between these
        2 cases because the second argument will be negative if we are nearer to
        the left endpoint.
      </p>
<p>
        Knowing this, we can rewrite our lambda expression to take advantage of this
        additional information:
      </p>
<pre class="programlisting"><span class="identifier">tanh_sinh</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span> <span class="identifier">integrator</span><span class="special">;</span>
<span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">,</span> <span class="keyword">double</span> <span class="identifier">xc</span><span class="special">)</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">x</span> <span class="special">&lt;=</span> <span class="number">0.5</span> <span class="special">?</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special">/</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="number">1</span> <span class="special">-</span> <span class="identifier">x</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">x</span> <span class="special">/</span> <span class="special">((</span><span class="identifier">x</span> <span class="special">+</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">xc</span><span class="special">)));</span> <span class="special">};</span>
<span class="keyword">double</span> <span class="identifier">Q</span> <span class="special">=</span> <span class="identifier">integrator</span><span class="special">.</span><span class="identifier">integrate</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">0.0</span><span class="special">,</span> <span class="number">1.0</span><span class="special">);</span>
</pre>
<p>
        Not only is this form accurate to full machine-precision, but it converges
        to the result faster as well.
      </p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
      Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
      Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde,
      Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle
      Walker and Xiaogang Zhang<p>
        Distributed under the Boost Software License, Version 1.0. (See accompanying
        file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
      </p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="de_tanh_sinh.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../double_exponential.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="de_sinh_sinh.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
